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Abstract. - We propose to detect the Mott insulator-superfiuid quantum phase transition in 
an array of coupled cavities by studying the polariton and photon fluctuations in a block of 
linear dimension M (in units of the lattice constant of the array). We explicitly show this for a 
one-dimensional array; the analysis can be however extended to higher dimensions. In the Mott 
phase polariton fluctuations are independent of the block size. In the superfluid phase they grow 
logarithmically with M, the prefactor being related to the compressibility of the system. In the 
case of photon fluctuations, the critical behaviour is encoded in the subleading scaling with the 
block dimension, while the leading behaviour is linear in M and non-critical. Our results have 
been obtained by means of the density matrix renormalization group numerical algorithm. 
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■ The recent suggestion [1-3] to realize Mott and super- 
fluid phases in arrays of coupled QED-cavities has stimu- 
lated a flurry of activity on strongly interacting photonic 
systems [4-12]. In the presence of randomness a glassy 
phase of polaritons was shown to appear in the phase dia- 
gram [5] . Coupled cavities can be engineered to behave as 
quantum simulators for a variety of interacting spin mod- 
els [7,9,10]; they can support soliton excitations [6,13]. 
A different realization of the strongly interacting regime 
with photons in an optical guide was proposed in [14], 
while polariton blockade effects were studied in a reso- 
nantly excited photonic quantum dot [15]. As candidates 
for simulating strongly interacting models, coupled cavi- 
ties present new characteristics as compared to other sys- 
tems, like optical lattices [16] or Josephson junction ar- 
rays [17]. Most notably, it is possible to access their local 
properties. 

In this paper we would like to exploit the local ad- 
dressability of cavity arrays to propose a method in or- 
der to probe the critical behaviour of the Mott insulator- 
superfiuid transition. Our suggestion is based on the study 
of the fluctuations in the number of photons and polari- 
tons. To our knowledge, this issue was not addressed so far 
in the literature. While it is well known how to distinguish 
a regime of photon blockade, how to detect the transition 
itself is, to a large extent, an unexplored problem. We 



study in details a one-dimensional array. Our conclusions, 
however, can be easily extended to higher dimensions. 

The Model - We suppose that inside each cavity a sin- 
gle two-level atom interacts with photons via a Jaynes- 
Cummings Hamiltonian [18]. While this situation is sim- 
pler to simulate, we do not expect, for the purposes of our 
work, any qualitative change if other atomic level struc- 
tures [1] are considered. The Hamiltonian for an array of 
L coupled cavities is then given by 



(1) 

where e denotes the transition energy between the two 
atomic levels, u> is the resonance frequency of the cavity, 
P is the atom-field coupling constant (e, u>, j3 > 0), and 
t is the inter-cavity photon hopping, that is assumed to 
be constant for all nearest neighbours and zero otherwise. 
The atomic and photonic raising/lowering operators are 
denoted as o~f, and {a\,ai} respectively, the subscript i 
indicating the lattice site. The total number of atomic 
plus photonic excitations (i.e., the number operator for 
polaritons) on the i-th site is given by ol = n? h + o-fa~ , 
where n^ h = a\ai is the photon number operator. Our 
analysis is restricted to the case of zero relative detuning 
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A = oj — e; we also work in the canonical ensemble with 
a fixed polariton density p = N/L = 1, N being the total 
number of polaritons in L cavities. 

The equilibrium phase diagram associated to the model 
in Eq. (1) is characterized by two distinct phases, with 
polariton Mott Insulating (MI) regions surrounded by the 
Supcrfluid (SF) phase. In the MI phase polaritons are lo- 
calized on each site due to the photon blockade [19], and 
there is a gap in the spectrum. A finite hopping rcnormal- 
izes this gap, which eventually vanishes at a critical value 
of t. For large hoppings excitations arc delocalizcd and 
the system enters the SF phase. 

We propose to detect the MI-SF transition by analyzing 
the fluctuations of the occupation number inside a block 
composed by a subset of M adjacent cavities. Since the 
number operator is the canonically conjugated variable 
with respect to the phase of the whole many-body wave- 
function, we expect that its fluctuations are strongly sup- 
pressed in the incoherent MI regime, and, by contrast, 
greatly enhanced in the coherent SF phase. The disper- 
sion of particle number on a given subsystem with M sites 
is quantified by the variance of the corresponding proba- 
bility distribution: 



Sn 2 a (M) 



igM 



(2) 



where a = pol/ph stands for polariton/photon fluctua- 
tions, and (•) denotes an average on the system ground 
state. This can be obtained from the two-point corre- 
lation functions of the related number operator: C" 



(i,j)£M ^ij 



{nfn°<) - (nf )(n?>, such that 8n 2 a (M) = £ 

In the present work we suppose that the system is in 
its ground state; we neglect decoherence, and assume that 
spontaneous photon emission and cavity loss characteris- 
tic times are much longer than the timescale over which 
the array can reach the ground state. We are aware that 
this is a strong assumption that may not be fulfilled in 
the experiments. We do not try to argue on this very 
important issue in the present work; rather we are inter- 
ested in describing a method that is suitable to detect 
the various (equilibrium and non-equilibrium) many-body 
phases of an array of cavities. We describe our proposal 
by using the (equilibrium) MI-SF transition; the essential 
features of the method can be as well applied to the non- 
equilibrium case. We will come back to this point in the 
concluding section. 

All the numerical data presented below have been ob- 
tained by means of the Density Matrix Rcnormalization 
Group (DMRG) algorithm with open boundary condi- 
tions. 

Polariton fluctuations - We first study fluctuations 
in the polariton number. The number of polaritons in- 
side a given subsystem can be experimentally measured 
by instantaneously switching off the effective polaritonic 
hopping (this can be achieved by changing the detuning 
A, such to increase the relative strength of the atom- 




Fig. 1: (Colour on-line) Variance 5ni y(M) of the probabil- 
ity distribution for the polaritonic occupation number inside 
a block of M cavities as a function of the inter-cavity photon 
hopping t, for fixed values of M. With our choice of parame- 
ters, the MI-SF transition point is at t* / j3 ~ 0.198 (red solid 
line). DMRG simulations have been performed for a system of 
L = 128 cavities, keeping a maximum of n max = 10 polaritons 
per site and a number m = 150 of states; blocks are formed 
starting from the center of the chain. 



field coupling with respect to the photon hopping). In 
this way, admitting radiative losses on long time scales, a 
quantum-jump picture immediately suggests that the po- 
lariton number is exactly given by the number of photons 
emitted from the selected region [20] . 

In Fig. 1 we show numerical data displaying 5rip ol (M) 
as a function of the photon hopping t in a system with 
L = 128 cavities; various curves correspond to different 
sizes M of the block inside the system. The MI-SF tran- 
sition point has been located numerically with high accu- 
racy at t*/j3 « 0.198 [5]. On-site fluctuations cannot re- 
veal critical features at the phase transition, because they 
correspond to a local property; whenever t ^ 0, the solid 
curve with circles corresponding to M = 1 has always non- 
zero values. Therefore, an unambiguous characterization 
of the phase boundary is impossible in this context. This 
was already pointed out in Ref. [21] for the on-site boson 
number fluctuations in the Bosc-Hubbard model, where it 
has been shown that the number probability distribution 
evolves from a Poissonian, in the noninteracting gas, to a 
sharply peaked distribution, in the insulator. Wc obtained 
very similar distribution probabilities for the on-site po- 
lariton number in our system, that are characterized by 
sub-Poissonian statistics. 

From Fig. 1 we notice that, by increasing M, in the MI 
phase fluctuations tend to be independent on the block 
size, while in the SF phase the dependence on M is evi- 
dent. In particular, we found that 8n^ ol (M) as a function 
of M saturates to a finite constant value in the MI phase, 
while it diverges logarithmically in the SF phase (see in- 
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set of Fig. 2, where we plotted the same curves of Fig. 1, 
but fixing t and varying M). The logarithm divergence is 
directly related to the fact that density correlation func- 
tions in the supcrfluid phase have a power-law decaying 
uniform term [22]: 



£rpol 



K 



(3) 



In the previous expression p is the particle (in this case 
the polariton) density and K is the so called Luttinger 
parameter [23], that is proportional to the square root 
of the compressibility of the system. Indeed, from the 
definition of dn^ ol (M), one immediately concludes that, 
in the SF phase 



8n 2 pol (M) 



In M . 



(4) 



At the SF-MI transition the coefficient K at integer den- 
sities is known to be equal to K c = 1/2 [23]. Therefore, 
for a polariton density p = 1, the logarithmic prefactor 
suddenly jumps from a value cq = 2/tt 2 just inside the SF 
phase, to zero in the MI, with a characteristic Kosterlitz- 
Thouless behaviour. 

A quantitative analysis of the crossover between the two 
different behaviours has been performed by fitting numer- 
ical data according to: 



«5^ ol (M)= C P ol ln 
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with A po1 and Cq°' as fitting parameters. The constant 
term A po1 is not important for our purposes, therefore we 
concentrate on the logarithmic term. In the main panel of 
Fig. 2 we plot the logarithmic prefactor c po1 as a function 
of t: the phase transition point in this situation can be 
quite clearly identified. In ID, in particular, a precise 
criterion for values of Cg 01 > 2/7T 2 in the superfluid phase, 
clearly identifies the transition point. 

Photon fluctuations - Wc now concentrate on the pho- 
ton number fluctuations, a quantity which is more directly 
measured in quantum optical experiments, where the sys- 
tem is typically in a non-equilibrium state between pho- 
ton/cavity decay and external pumping. The state of the 
polariton field is usually retrieved by detecting and char- 
acterizing the light emission [24,25]. 

In this case the situation is quite different as compared 
to the case of polariton fluctuations. Even at zero hop- 
ping the on-site photon number is fluctuating. For ex- 
ample, in the case of zero relative detuning (A = 0) 
and deep in the MI regime, we have (n ph ) w 1/2 and 
(K h ) 2 ) « 1/2- Therefore, even for a perfect insulator, 
the variance Sn 2 h (M) of the photon number distribution 
inside a block of length M is proportional to the block 
size: 5np h (M) w M/4. The onset of the superfluid be- 
haviour can then be sought in the deviations from this 
linear growth. These deviations are due to the raising of 




Fig. 2: (Colour on-line) Best fitting parameter Cq° of Eq. (5) 
for the variance 8n^ ol (M) as a function of the hopping t. The 
horizontal dashed line indicates the estimate of the logarithmic 
prefactor at the transition point: Co ~ 2/tv 2 . In the inset we 
plot 8np 0l as a function of M (symbols) for some fixed values 
of t: from bottom to top t = 0.12, 0.16, 0.2, 0.24, 0.28, 0.32; 
continuous lines are logarithmic fits of the corresponding data. 
In numerical fits we dropped the first and the latter 10 points 
(i.e., we kept values for M € [11, L/2 - 10]). 



correlations between distant sites on increasing the hop- 
ping strength. Therefore, in the scaling ansatz for the 
photon fluctuations one should also include a term that is 
linear in the block dimension, i.e. 



Snl h (M) = aM + cf In 



L . 




— sin 




7T 





+ A ph . (6) 



Let us analyze the properties of the first two terms in the 
r.h.s. of Eq. (6). The behaviour of a as a function of 
the hopping is shown in Fig. 3. Starting from the value 
a = l/4att = 0, the coefficient of the linear term de- 
creases on increasing the hopping. This behaviour signals 
the fact that correlations between photon fluctuations at 
different locations start to develop. Although at t ~ 0.2 
there is an indication of the change in the ground state 
properties, the coefficient a seems not to be an appropri- 
ate indicator of the location of the critical point. As for 
the case of polariton fluctuations, also in the case of pho- 
tons one should look at the coefficient of the logarithmic 
term. This is shown in Fig. 4: here the critical behaviour 
is easily identified. The value of the fluctuation at the 
transition is not universal as in the polariton case; the 
reason is that photons are not the genuine excitations of 
the system, therefore the considerations leading to Eq. (4) 
with K c = 1/2 do not apply. The ratio between the co- 
efficients of polariton and photon fluctuations is shown in 
the inset of Fig. 4. A marked change at the transition can 
be noticed. In our opinion, however, this is not obviously 
related to the critical point; rather we expect that the ra- 
tio is a non-universal feature depending on the details of 
the model and on the value of the couplings. 
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0.26 




Fig. 3: (Colour on-line) Coefficient of the linear term in the 
photon fluctuations, see Eq. (6), as a function of the photon 
hopping t between the cavities. As it is evident from the plot, 
a carries no information about the critical behaviour (as a 
reference, we indicated with a red bar the location of the critical 
point). 



0.5 
.ph 



0.4 



0.3 



0.2 



0.1 



2.6 

2.4 

0-02.2 
O 

2 

"3 

■JO 1.8 



V 



0.1 0.2 



0.3 ,o 0.4 

r 



0.1 



0.2 



0.3 



0.4 



Fig. 4: (Colour on-line) Prefactor Cq 11 of the logarithmic term 
in Eq. (6) as a function of the hopping. The behaviour is 
very similar to the corresponding prefactor for the polariton 
fluctuations. The ratio between Cq o1 and Cq , shown in the 
inset, is non-universal. 



Indication of the critical behaviour from photon correla- 
tions may also be obtained without resorting to a scaling 
analysis as a function of the block size M. This can be 
achieved by studying the second derivative of photon num- 
ber fluctuations with respect to the block size, evaluated 
for blocks of half the length of the system size, that is 
d\j [<5rip h (M)] \ M _ L / 2 ' Although this approach might be 
slightly less accurate, it can give an interesting insight into 
the critical region. 



We first consider d\j [Sn^ h 



(M)] \ m _t /, as a function of 



\M=L/2 

the hopping parameter t. Numerical data in Fig. 5 show 
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Fig. 5: (Colour on-line) Second derivative of photon number 
fluctuations with respect to the block length, evaluated at half 
the length of the system size, as a function of the hopping t. 
Inset: logarithmic prefactor of polariton number fluctuations, 
as a function of the second derivative of photon number fluc- 
tuations (data are the same of the main panel and of Fig. 2). 



that, like for the polaritonic number fluctuations, also this 
quantity can be used to characterize the MI-SF transition 
for finite sizes. In addition, we can identify a behaviour 
that is very similar to the one of the logarithmic prefac- 
tor Cq°' of Eq. (5) for the polariton number fluctuations. 
Namely, the correlation between these two quantities is 
nearly perfect: apart from a proportionality factor, which 
depends on the system size, the two curves in the main 
panels of Figs. 2-5 are exactly the same. This is shown in 
the inset of Fig. 5 where, for a given value of the photon 
hopping t, we plot the corresponding values of Cq o1 and of 
d\j [<$np h (M)] \ M=L i 2 in the two axes, thus displaying a 
perfect linear correlation. 

The dependence of the photon fluctuations on the size 
is however more complicated than that of Srip ol . In Fig. 6 



we plot d 2 M [<5rip h (M)] | M=i / 2 (its absolute value, actually, 
since it is always negative) as a function of the system size 
L, for several fixed photon hoppings t. While the logarith- 
mic prefactors Cg 01 and Cg h are independent of the system 
size, this is not the case for the second derivative of photon 
number fluctuations, which asymptotically drops to zero 
as a power-law with L. More specifically, in the free pho- 
ton limit (t — > +oo) it decays as L~ 1 , while for finite values 
of the photon hopping and a sufficiently large system size, 
there is a crossover to a L~ 2 behaviour. This follows from 
two qualitatively different decays of the photon number 
correlation function C 1 ?' 1 = (n^ r n^ h ) — (n?^ r )(nf h ): apart 
from open-boundary effects, in the first regime C^' 1 ps 1/L, 
while in the second case C? h ps 1/r 2 . We carried out a 
numerical analysis of the photon correlations C^ h in our 
QED-cavity array system, and explicitly found these two 
distinct behaviours. 

It has also been possible to derive an analytic estimate 
of the crossover scale by exploiting a mapping of the cavity 
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Fig. 6: (Colour on-line) Second derivative of photon number 
fluctuations with respect to the block size M evaluated at 
M = L/2, as a function of the system size L. The two straight 
blue lines denote power-law behaviours L~ 2 . DMRG pa- 
rameters are the same as in Figs. 1 and 2. 



array with a Bose-Hubbard (BH) model, and studying its 
corresponding boson correlators C™. It should be kept 

in mind that the formal analogy between Cf^ and Gy H 
has to be considered only at a qualitative level, since the 
equivalence of the two models becomes exact only in the 
limit of a large number of atoms per cavity [1]. 
The BH Hamiltonian is defined by 



(7) 



where {bj,bj} are the boson creation/annihilation oper- 
ators, rii = b\bi the boson number operator, and with 
C BH wc indicate the corresponding correlation functions. 
When the depletion of the condensate is not too great 
(J 3> U), we can employ the Bogoliubov approximation, 
which consists in replacing the boson creation and anni- 
hilation operators at a given site j by a c-number z% E C 
plus a fluctuation operator (3j . Only terms at most of the 
second order in (3j are considered when diagonalizing the 
Hamiltonian [26-28]. Within this framework and for pe- 
riodic boundary conditions, we locate the crossover scale 
between the L _1 and L~ 2 behaviour by 



U/J 



< 



2tt 2 
n Q L 2 



(8) 



with no being the density of the boson condensate. When- 
ever the inequality is valid, the spectrum of the Bogoliubov 
quasi-particles at small k is quadratic, and correlations be- 
have as for free bosons: Cj 311 «1/X. Conversely, if it does 
not hold, the spectrum becomes linear at small k, and the 
free-boson limit fails; in this regime we recover the decay 
C ? BH ~ l/r 2 for large r, typical of the polariton number 
correlations [22]. 
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Fig. 7: (Colour on-line) Number fluctuations of photons C£ h 
in an array of cavities (black curves - DMRG data) and of 
bosons C^ H in BH model (red curves - Bogoliubov ansatz); 
dashed blue curves display a behaviour C r ~ r~ 2 . The two 
panels correspond to different values of the hopping, that is 
assumed to be equal in the two models t = J. In the insets 
we plot the deviation a of correlations for the two models, as 
a function of U ; the red curve in each panel shows the best 
fitting curve of DMRG data, corresponding to the value of U 
which minimizes a. We analyzed chains of L — 128 sites, and 
evaluated correlations starting from the central site. 



A proper quantitative comparison between photon num- 
ber correlations in the coupled cavity system and boson 
number correlations in the BH has been obtained by nu- 
merically solving the BH Bogoliubov equations for a sys- 
tem with open boundaries, and it is shown in Fig. 7. The 
hopping is chosen to be equal in the two models (t = J); 
the polariton/boson density is p = 1. DMRG data are 
fitted by varying the on-site repulsion U in the BH and 
by minimizing the squared differences between the two 
curves: a 2 (U) = J2 r \ C ? h - C™\ 2 . Of course, if t/0 is 
too small the Bogoliubov approximation fails, and we do 
not find an appropriate minimum. We remark that corre- 
lations are always negative, thus indicating photon anti- 
bunching, as it is revealed by on-site number distributions 
with negative values of the Mandel parameter [18]. 

Conclusions - In this paper we studied photon and po- 
lariton fluctuations in coupled cavities. By a suitable scal- 
ing of the fluctuation detected over a region of the sample, 
we showed that it is possible to extract the critical prop- 
erties of the superfluid to Mott insulator phase transition. 
The analysis presented here for one-dimensional systems 
can be extended to higher dimensional lattices. We remark 
that we confined our study to the equilibrium case; an in- 
teresting point that remains to be addressed is to consider 
the effects of photon and cavity losses, which may lead 
to new non-equilibrium effects in the physics of such sys- 
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terns. We do however expect that the method proposed 
here will hold also in the non-equilibrium case. The reason 
is that our approach ultimately distinguishes the system 
being incompressible or not. While the general princi- 
ple also holds for driven systems, the detailed behaviour 
leading to the logarithmic scaling may be modified. As a 
consequence, the most appropriate scaling ansatz for the 
non-equilibrium cases has to be checked. 

* * * 
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